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O ABSTRACT 

a: 

The maps presented in Paper I are here used to infer the variation of the column 
\ densities of HCO + , DC0 + , N2H + , and N2D + as a function of distance from the 

^ ' dust peak. These results are interpreted with the aid of a crude chemical model 

C$ • which predicts the abundances of these species as a function of radius in a spherically 

symmetric model with radial density distribution inferred from the observations of 
dust emission at millimeter wavelengths and dust absorption in the infrared. Our main 
observational finding is that the A r (N2D + )/A r (N2H + ) column density ratio is of order 
0.2 towards the L1544 dust peak as compared to N(DCO + )/N(UCO + ) = 0.04. We 
conclude that this result as well as the general finding that N2H + and N2D + correlate 
well with the dust is caused by CO being depleted to a much higher degree than 
molecular nitrogen in the high density core of L1544. Depletion also favors deuterium 
enhancement and thus N2D + , which traces the dense and highly CO-depleted core 
nucleus, is much more enhanced than DCO + . Our models do not uniquely define 
the chemistry in the high density depleted nucleus of LI 544 but they do suggest 
that the ionization degree is a few times 10 -9 and that the ambipolar diffusion time 
scale is locally similar to the free fall time. It seems likely that the lower limit which 
one obtains to ionization degree by summing all observable molecular ions is not a 
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great underestimate of the true ionization degree. We predict that atomic oxygen is 
abundant in the dense core and, if so, H30 + may be the main ion in the central highly 
depleted region of the core. 



Subject headings: ISM: individual (L1544) - ISM: dust, extinction - ISM: molecules 



1. Introduction 

The ionization degree x(e) (= n(e)/n(H2), with n(e) and n(H2) the electron and H2 number 
density, respectively) in dense molecular clouds plays a key role in determining the initial conditions 
which precede the collapse to form a star (e.g. Shu, Adams & Lizano 1987). If the magnetic field 
threading the dense gas is sufficiently large to prevent immediate collapse, ambipolar diffusion 
of neutrals across field lines can lead to a situation where a dense core of gas is gravitationally 
unstable. The time scale for such ambipolar diffusion is proportional to the ionization degree and 
it therefore becomes of interest to develop methods of estimating the ionization degree based upon 
the abundances of various species which trace the ionization in the dense gas. 



In an earlier article ( |Casclli et al. 1998b , hereafter CWT98), we studied the possibility of 



using the [DCO + ]/[HCO + ] abundance ratio as a tracer of the ionization degree in dense molecular 



gas, following previous work by Wootten, Snell & Glassgold 1979 and Guelin, Langer & Wilson 



1982 , We used the time-dependent models of Lee et al. 1996| to study the expected behavior 



of a variety of species thought to be relevant in this context and compared with the abundances 



measured by Butner et al. 1995. An analogous study carried out by Williams et al. 1998 



reached very similar results. We confirmed earlier work ( Dalgarno fc Lepp 1984j ) which showed 



that a critical parameter in addition to the ionization degree which determines the fractionation 
of deuterated species (and hence abundance ratios such as [DCO + ]/[HCO + ] ) is the degree of 
depletion of carbon and oxygen bearing species onto the surfaces of interstellar dust grains. 
Measuring the ionization degree in reliable fashion therefore requires the use of observations to 
solve for both x(e) and depletion. One possible approach to this problem discussed by CWT98 



was to use the cyanoacetylene (HC3N) abundance as a "depletion indicator" (see also Ruffle et al. 



1997[) . A more detailed analysis on the possibility to use HC3N to estimate the degree of depletion 
in a dense core is currently under way (Comito et al., in prep.). 

The CWT98 and |Williams et al. 1998] studies used data with angular resolution of the order 
of an arc minute and did not search for spatial variations in the ionization degree and hence in 
the [DCO + ]/[HCO + ] ratio. Studies of variations in the deuterium fractionation and in the degree 
of ionization across cores have been previously done in active star forming regions such as the 



R Coronae Australis molecular cloud ( Anderson et al. 1999| ), and in massive cores located in 



regions that are currently forming stellar clusters ( [Bergin et al. 1999| ). In the former case, the 
electron fraction was found to decrease near the cluster of newly born stars, whereas in the latter 
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no differences in the ionization fraction were detected between cores with and without associated 
stars. No changes in the ionization fraction was found either between the edge and the center of 
the mapped regions. However, no similar studies have been performed so far in the direction of 
pre-star-forming cores in relatively quiescent zones of our Galaxy, which give the opportunity 
to investigate the initial conditions in the star formation process presumably not triggered by 
external agents. 

In any model of a pre-protostellar cores, one expects to observe a concentration of material 
around the region of highest density. As the density becomes higher, one expects the degree 
of ionization to decrease, the degree of depletion onto dust grain surfaces to increase, and, 



consequently, the fractionation of deuterated species to be enhanced (e.g. Roberts Millar 



2000). In fact, one might naively expect that deuterated species correlated with high depletion 
might trace the phase prior to the formation of a protostar. With this in mind, we undertook 
observations using the IRAM 30-m telescope of a variety of species towards the dense core L1544 



in the Taurus molecular cloud which is thought (e.g. Tafalla et al. 1995 , Ciolck and Basu 2000| ) 



to be in a phase shortly prior to that at which collapse occurs. The data and a kinematical 
study of this region are presented in Casern et al. (2001, hereafter Paper I). Here we present 
our determination of the electron fraction across L1544 and show that x(e) at the center of the 
core - where the H2 number density is ~ 10 6 cm -3 - is about 2xl0 -9 , about five times smaller 
than expected if the electron fraction is due to cosmic-ray ionization alone and the cosmic-ray 
ionization rate ( has its "standard" value (x(e) = 1.3x 10" 5 n(H 2 )~ a5 ; cf. pVlcKee 198S|) . 



In section 2, we show our determination of physical and chemical parameters of importance 
for electron fraction estimates. The determination of electron fraction is presented in Section 3. 
Finally in Section 4, we discuss the consequences of our measurements for the understanding of 
high density cores like L1544. 



2. Determination of physical and chemical parameters 

An accurate estimate of electron fraction implies measurements of (i) the column density of 
the observed molecular ions, (ii) the temperature as well as (hi) the density structure of the core, 
and (iv) the amount of CO depletion. In the following sections we will describe the above points 
separately. 



2.1. The temperature across the core 

Because of its low dipole moment, CO is likely thermalized at the typical core densities and 
it can be used to measure kinetic temperatures. |Tafalla et al. 1998 used the J = 1— >0 transition 
of 12 CO, 13 CO and C 18 to determine the excitation temperature of the CO emitting gas across 
L1544 and found T ex ~ 12 K, in good agreement with previous kinetic temperature estimates 
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based on CO and NH3 observations (~ 10 K; Myers et al. 1983| , [Ungerechts et al. 1982| , Benson 
fc Myers 1989 ). We repeated a similar analysis using C 18 O(l-0) and C 17 O(l-0) data and found 
an average value of T ex = 10.1±0.1 K . Of course, CO is not a good tracer of gas temperature 
at densities ^ 10 5 cm -3 , where it is strongly depleted (CWT99; see Sect. 2J3). A recent Monte 
Carlo radiative transfer analysis of the two inversion transitions of NH3, observed at the Effelsberg 
antenna (Tafalla et al., in prep.; 40" half power beamwidth), has also shown that Tj< ~ 10 K with 
no significant variation across L1544. We assume on this basis a constant temperature of 10 K in 
the rest of our analysis. However, model calculations of the dust temperature in pre-stellar cores 
( |Zucconi et al. 2001 , Evans et al. 2001 ) predict a fall off in dust temperature from ~ 14 K at 
the edge to ~ 7 K at the center. This is likely to be reflected also in the behaviour of the gas 
temperature at least for densities above 10 5 cm" 3 . Although for our purposes we can neglect such 
findings (our simple chemical model presented in Sect. 3.2 is not significantly affected by a change 
in temperature from 10 to 7 K), the predicted temperature gradient affects our estimates of the 
volume density from N2H + , as discussed below. 



2.2. The density across the core 



From our molecular line data we can deduce good number density estimates only at the 
L1544 "molecular peak" and in a few adjacent positions (see Sect. 2.2.1, p.2.2 ). Therefore, for our 



electron fraction estimates across L1544 (Sect. 3.2) we used the density profiles obtained from 
1.3mm continuum observations ( [Ward-Thompson et al. 1999| , hereafter WMA) and ISOCAM 
absorption in the mid-infrared ( Bacmann et al. 2000| , hereafter BAP). These density profiles are 



described in Section 3.2.1. 



A number of line intensity ratios measured by us are sensitive to density. Since the line ratios 
depend on collisional rates, one in principle derives n(H2) independently of assumptions about 
abundance. However, in a situation where the density varies rapidly along the line of sight (as 
suggested by the dust emission and absorption), the density which one derives will be a weighted 
mean. In the following, we describe some estimates we have made using different tracers. 



2.2.1. Density inferred from DCO+ (3-2)/DCO + (2-1) 



Our observations of DCO + (3-2) and (2-1) allow estimates to be made of the hydrogen 
number density in the regions responsible for the emission of these lines as well as estimates of the 
DCO + column density. We have done this using an LVG statistical equilibrium code similar to 
that described e.g. in Walmsley (1987). The approach is essentially identical to that discussed by 
Butner et al. 1995} The collisional rates were taken from those calculated for collisions of para-H2 
with HCO + by Monteiro (1985) assuming a temperature of 10 K. At this temperature, the rates 
are essentially identical to the more recent rates of Flower 199S . We assumed a dipole moment of 
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3.9 Debye following Botschwina 1989| (see also Williams et al. 1998 ). 



In Fig. 0, we show the results obtained in this manner superimposed upon the half-maximum 
contours of the DCO + (2-l) and (3-2) maps (see Paper I for details on these data) and the 1.3mm 
continuum map of WMA. We see a definite tendency for the highest densities derived from DCO + 
to be close to the dust emission maximum. We see also that the values found for n(H2) towards 
the dust peak are close to those derived by WMA from their observations (~ 10 6 cm -3 ). However, 
there are a number of caveats of which one should be aware. 

One of these is that the DCO + (2-l) line is almost certainly optically thick towards the core of 
L1544. This is suggested by both our LVG results themselves and the comparison at offset (20,-20) 
with D 13 CO + (see Paper I) which suggests an optical depth of order 4. Our density estimates are 
thus sensitive to the geometry and velocity gradients within L1544. We note additionally that 
our DCO + (3-2) map cannot be convolved to the 2-1 resolution (1.5 times worse) because it is 
undersampled (grid of 1.7 times HPBW). Nevertheless, our results suggest that we are sampling 
in DCO + the high density inner core of L1544. 



2.2.2. Density derived from N2H + and N2D + 

Our N2H + and N2D + observations allow an independent estimate of the density. This is 
of special interest because, as we have seen in Paper I, the N2H + and, in particular, the N2D + 



integrated intensity maps appear to trace best the dust emission (see also Bergin et al. 2001). 
Another characteristic of these species is that, due to the hyperfine splitting, we can determine 
directly the optical depth in several transitions utilising the relative intensities of the hyperfine 
satellites. This reduces considerably the errors in our determinations of level column densities 
for, say, N2H + , as compared to DCO + . Finally, it is of some importance that we have detected 
N2H + (3-2) and N2D + (3-2). This is because we expect such transitions to be essentially collision 
dominated in the sense that the integrated line intensity has, for a given temperature, a one-to-one 
correspondence with the integral of the product of the H2 density and the N2H + density along 
the line of sight. We call this integral / n(H2)n(N2H + )ds the N2H + collision measure CM(N2H + ) 
in the following. Essentially, the (3-2) line intensities are proportional to CM(N2H + ) in the 
limit where collisional deexcitation is negligible with respect to radiative transitions which for 
N2H + (3-2) at 10 K for example implies n(H2) less than 10 6 cm -3 . 

We have developed an LVG code for N2H + (N2D + ) analogous to that for DCO + . We use the 
collisional rates of Green(1975) for collisions of N2H + with He and multiply them by 1.4 to take 
account of the differing mass of H2. We also have adjusted the escape probability to account for 
hyperfine splitting in the 1-0, 2-1, and 3-2 transitions of N2H + . We first use our determination 
of the optical depth in the N2H + (l-0) line at offset (20, -20) to obtain an estimate of the N2H + 
column density and find iV(N 2 H+) = 1.5 10 13 cm" 2 at offset (20, -20). Then, we determine 
CM(N2H + ) from N2H + (3-2) data, which have been averaged in a 3x3 grid of positions spaced 
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by 10" and centered at (20, -20) (see Tab. 3 of Paper I, second row) in order to simulate a 20" 
beam, similar to the 1-0 observations. From the CM(N2H + )/iV(N2H + ) ratio we find a mean 
density of 1.3 xlO 5 cm -3 , smaller than the value derived above from our DCO + measurements 
and also smaller than the central densities suggested by the results of WMA and of BAP. These 
discrepancies may be partly due to the assumption of a constant temperature of 10 K, in particular 
in the dense highly CO-depleted nucleus of L1544, best traced by N-bearing molecules, where 



Zucconi et al. 2001 derive T ~ 7 K 



We can in analogous fashion use our N2D + observations to obtain a density estimate. At the 
center of L1544 assuming a temperature of 10 K, we find a mean density CM(N2D + )/./V(N2D + ) 
= 1.510 18 /4.910 12 = 310 5 cm . This is a factor of 2 higher than the N2H + estimate above 
which might perhaps reflect the fact that the N2D + emission is more compact than that of N2H + . 
However, it is still a factor of 3 smaller than the density deduced from dust continuum and DCO + 
observations, again suggesting a possible temperature fall off in the core center. 



2.3. Depletion of CO 

The depletion factor /d is defined by the ratio between the "canonical" fractional abundance 
of CO (9.5xl0^ 5 , [Frerking et al. 1982| H) and the observed x(CO) = JV(CO)/N(H 2 ). To determine 



/d at the dust peak position, CWT99 compared the C O integrated intensity with the 1.3mm 
continuum dust emission flux density from WMA. They found that /d ~ 10 at the dust peak, 
and that observations are consistent with a model where CO is condensed out onto dust grains at 
densities ~ 10 5 cm -3 . The corresponding radius of the region where CO is severely depleted is ~ 
6500 AU and the depletion causes 2.3 M of gas to be lost to view in molecular line emission. 

Here we extend this analysis to all the other positions where C 17 has been observed. To 
do this we divided the 1.3mm map by the CO (1-0) map. This division implies three steps: (i) 
convolve the 1.3mm map with a 22" beam (the HPBW of the C 17 observations) and reproject 
it to have the same coordinates as the C 17 map; (ii) make a regularly sampled C 17 map; (iii) 
divide one map by the other. The result of this operation can be translated into a depletion 
factor map. We somewhat arbitrarily assume "normalcy" to be given by [C 17 0]/[H2] = 4.8xl0 -8 
(Frerking et al. 1982). Assuming an excitation temperature for C 17 of 10 K and a 1.3 mm dust 



grain opacity k u = 0.005 cm g (WMA) , we find that the depletion factor /d can be expressed 



as: 



_ 5(1.3 m n, mJy /22») 
W C 17 ( KkmS " > 



We note that the only direct measurements of the CO abundance give larger values than those reported by 
Frerking et al. 1982) by a factor of about 5 (Lacy et al. 1994), implying larger depletion factors than found in this 



work. 
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Here, ^q17q 1S the integrated intensity ( over all hyperfine components) of C 17 O(l-0) and 
5(1. 3mm, mJy/22") is the flux density at 1.3mm in a 22" beam. In those positions where only 
C 18 was observed, we used / D = 0.085 5(1.3mm, mJy/22")/T^ c i 8o (Kkms" 1 ). 

In Fig. [2] the depletion factor map is shown together with the 1.3mm continuum map from 
WMA. There is a fairly good correspondence between the /d and 1.3mm contours as one might 
expect given the "flat" nature of the C O and C 18 maps. We note however that the peaks are 
not coincident and there is a projected distance between the /d and the dust peaks is 13" (about 
half a beam). This appears significant but one must remember that the C 18 and C O maps are 
probably dominated by emission from gas of density a few times 10 4 cm -3 to the foreground or 
background of the dust emission core. Thus structures unrelated to the pre-protostellar core may 
influence our maps. 



2.4. Molecular ion column densities 

The methods used to estimate column densities of the observed species, which give same 
results within a factor of 2, are reported in appendix [A|. Here we show the results. 



2.1.1. DCO+ and HCO+ 



The statistical equilibrium calculations used to estimate the H2 density (see Sect. ^2) can 
also be applied to a determination of the column densities of DCO + and HCO + . As in the case of 
density determination, our errors grow rapidly if the observed lines are optically thick. Towards 
the peak of L1544, we have used the optically thin species HC 18 + and D 13 CO + to determine 
HCO + and DCO + column densities and C 17 to determine the CO column density. We have 
considerable confidence that in these three cases, the observed lines are in fact optically thin (see 
discussion in Sect. A.l). For C 17 0, we can check (see e.g. CWT99 ) optical depth using the 



relative intensities of the hyperfine satellites. 

We nevertheless for most of the following discussion use H 13 CO + and DCO + column densities 
derived from transitions which are optically thick towards the dust emission peak and in the 
immediate vicinity. The reason for this is the greater extent of our maps in these lines. The 
justification is that our checks with the optically thin variants yield results which are consistent 
with those presented here. Thus we have checked our H 13 CO + column densities using HC 18 + 
where available and our DCO + column density using the D 13 CO + measurement at (20,-20). We 
conclude that our column density estimates are accurate to within a factor of 2 also when the 
optical depth is high (see appendix [A. 1 for details). 



Bearing this in mind, we present in Fig. |3| our DCO + and HCO + column density maps. One 
sees in the first place that these differ considerably from the integrated intensity maps. While 
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for the reasons discussed above, the column density maps (or ratio maps) should be treated with 
caution, we consider them a better estimate than simply scaling the integrated intensity maps. 
We also show in Fig. [3| the inferred ratio of column densities iV(DCO + )/iV(HCO + ). DCO + is 
clearly more concentrated towards the structure seen in dust emission by WMA than is H 13 CO + . 
Considering the column density ratio map, we see that [DCO + ]/[HCO + ] appears to vary by almost 
an order of magnitude between the centre and the edge of the map (roughly an arc minute or 0.04 
parsec). Away from the peak, the inferred iV(DCO + )/-/V(HCO + ) ratios are considerably lower 
than found in the survey of Butner et al. 1995| . This is somewhat surprising because we expect to 
be observing higher density gas on average with greater degrees of depletion than in the arc minute 
resolution [Butner et al. 1995| survey. The |Butner et al. 1995| estimates for jV(DCO + )/jV(HCQ + ) 
vary between 0.02 and 0.07 whereas our results for L1544 vary from 0.04 near the dust peak to 
values a factor of 10 lower at the edge of our map. 



24.2. N 2 D+ and N 2 H+ 

The N 2 D + and N 2 H + column density maps, together with A^(N 2 D + )/./V(N 2 H + ), are also 
shown in Fig. ||. Our N 2 D + observations as noted previously (see also Paper I) show that N 2 D + 
traces the dust peak emission better than any other species. It also turns out that the values we 
infer for iV(N 2 D+)/iV(N 2 H+) are much higher than for iV(DCO+)/iV(HCO + ). At offset (20,-20), 
we infer iV(N 2 D + )/./V(N 2 H + ) = 0.24±0.02. This is a factor of 4 higher than our maximum value 
for A r (DCO + )/A r (HCO + ). We conclude from this that the degree of deuterium enhancement 
increases towards the dust emission peak of L1544. 



3. The electron fraction x(e) 

CWT98 used a simple chemical model to find analytic expressions which allow one to directly 
estimate x(e) (and the cosmic ray ionization rate £) once 7V(DCO + )/iV(HCO + ), /d, and n(H 2 ) 
are known from observations (see their eqns. 3 and 4). However, these expressions furnish upper 
limits of x(e) and £ (see Section |3~T| ) . In fact, in the simple chemical scheme used to derive CWT98 
analytic expressions, dissociative recombination of molecular ions onto negatively charged dust 
grains and gaseous reactions involving atomic deuterium (see e.g. Dalgarno fe Lepp 1984| ) are not 



taken into account. Moreover, this chemical model is only valid for homogeneous clouds, whereas 
LI 544 clearly has a density structure and a varying amount of depletion and molecular abundances 



along the line of sight (see Sect. 2A). In fact, CWT98 eqns. (3) and (4) have been applied to data 
with a significantly lower (factor of ~ 6) spatial resolution compared to the present observations, 
so that any gradient was smoothed out in the beam and resultant A r (DCO + )/iV(HCO + ), /d, and 
n(H 2 ) should be considered as "average" quantities for the cores. The same problem is met if 
"pseudo-time dependent" chemical codes are used. The density structure as well as differential 
depletion of molecular species has to be taken into account to reach a satisfactory agreement with 
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observations (see also Aikawa et al. 2001). 

In this section, we discuss what one can say about the ionization fraction in L1544. One may 
obtain a lower limit to the ionization fraction x(e) summing the abundances of observed molecular 
ions in L1544. This limit may be close to the real ionization fraction, if refractory metals of low 



ionization potential such as Mg and Fe are significantly depleted (see e.g. Fig. 2 of Caselli et al 
l998bD. 



We first briefly discuss the limits one can place on ionization degree based upon the observed 
molecular ion column densities. 



3.1. Lower limits for x(e) 

Lower limits of the electron fraction across L1544 can be obtained by simply summing up 
column densities of the observed molecular ions: 

^ iV(HCO+)+iV(DCO + ) + 7V(N 2 H+) + 7V(N 2 D+) iV(HCO+) 

x{e)l - mh) WZT' (2) 

x{e)i has been obtained for each observed position toward L1544 based upon the column density 
determinations from Sect. [2.4| and the column densities of WMA. We find in this way a lower 
limit at essentially all positions where H 13 CO + was observed of : < x(e)i > ~ lxl0~ 9 with little 
variation. 

These estimates do not however really sample the high density nucleus of LI 544 and in order 
to do that, we use the model discussed in the following section. 



3.2. x(e) from chemical models 

Deriving the behavior of the electron fraction as a function of radius is complicated however 
by the strong abundance gradients due to depletion. These cause observed quantities such as e.g. 
the HCO + column density to be strongly influenced by background and foreground emission in 
the observed lines. This is the main rationale for the simple model introduced below. 



3.2.1. Model description 

The results of WMA and BAP show rather clearly that the L1544 "core" itself contains a high 
density "core" where the density reaches values as large as 5 10 5 cm -3 over a volume of radius 
several thousand AU. We have shown previously (CWT99) that CO is depleted by at least an 
order of magnitude in the central high density peak. In this section, we consider the behavior of 
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the ions and of the electron density in the depleted region. We have done this with an extension 
of the model briefly described by CWT99 and used to fit the CO depletion. 

We assume for this purpose a spherically symmetric density distribution with a constant 
molecular hydrogen density of (i) 5.4 x 10 5 cm -3 out to a radius rg at = 2900 AU followed by a 
1/r 2 fall off out to a cut-off radius of 10000 AU, which simulates the results of BAP based on the 
extinction measured with ISOCAM against the emission of small particles at a wavelength of 7 
fj,m; (ii) 1.3xl0 6 cm -3 out to a radius rfj at = 2500 AU followed by a 1/r 2 fall off out to a cut-off 
radius of 13000 AU, as found by WMA on the basis of 1.3mm observations of dust emission. 
The H2 column and volume density inside rji a t determined by BAP are a factor of 3.8 and 2.4, 
respectively, lower than those found by WMA (see however the comments in Table 2 of BAP). 
It is clear that this approximation is somewhat crude and indeed that the distribution departs 
considerably from spherical symmetry. Moreover, the recent study of [Evans et al. 2001 has shown 



that the temperature drop in the L1544 core nucleus to about 7 K (in agreement with ^ucconi"et 



al. 2001|) impl les a steeper density profile than in the isothermal case of WMA and BAP. However, 



we believe that our simple model suffices to explain many of the observations. 

As in CWT99, we have followed the time-dependent depletion of CO neglecting dynamic 
effects. That is to say, we assume the depletion time-scale to be more rapid than dynamical 
time-scales and assume initially a "canonical" CO abundance of 9.5 x 10~ 5 relative to H2. 
However, we also consider in this study the recycling of CO back into the gas-phase using the 
cosmic ray desorption rate proposed by Hasegawa fc Herbst 1993] and a CO adsorption energy 



Ed(CO) to the grain surface of 1210 K ( Hasegawa et al. 1992; ). This is for a surface of Si02 



while for binding with a CO mantle, Sandford and Allamandola (1990) find £x>(CO) to be 960 K. 
This is critical for the model because we have, following [Hasegawa Sz Herbst 1993| , assumed a 
time-scale for cosmic ray desorption t cr of species X given by : 

l/t cr {X) = 310- 7 Ci7 x exp[-E D (X)/70] (3) 

Here, £17 is the cosmic ray ionization rate in units of 10~ 17 s _1 and Ed(X) is the adsorption energy 
in K of species X. Rough estimates suggest that "chemical" desorption due to molecular hydrogen 
formation ( Willacy &: Millar 1998| ) appears to be less effective. Indeed, a recent theoretical work 



by Takahashi fc Williams 200C| has shown that the chemical desorption of CO can occur on 



small grains with size less than 20 A, but it is negligible on larger grains. A crucial point in 
our discussion moreover is that we also consider the analogous depletion of molecular nitrogen 
assuming an initial N 2 fraction of 7.5xl0~ 5 flMeyer et al. 19971 ). This is a maximal value in that 



it is based on nitrogen depletion measurements in diffuse clouds. Critical for our model is also 
that we adopt a N2 adsorption energy of < 787 K (= 0.65xI?d(CO); [Bcrgin fc Langer 1997 based 



upon calculations of [Sadlej et al. 1995 ) considerably lower than that for CO. The exponential 



dependence of the cosmic ray desorption rate above causes more effective depletion of CO and 
gives rise to a layer where N2 is the most abundant gas phase species containing heavy elements 
and where N2H + is an abundant ion. 
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In a model of this type, one can write for the abundance n(X) of species X (where X can be 
CO or N2) relative to H2 that : 

n{X) = n(X,oo) + (n tot (X) - n(X, 00)) exp-(t/t ). (4) 

Here, ntot(X) is the total abundance of X in both gas and solid phase , n(X, 00) is the steady 
state gas phase abundance given by ntot(X)/(l + t cr /td ep ) and tdep is the depletion time scale 
1/(S n gT a gr v (X)) (where S is the sticking coefficient, n gr is the grain number density, v(X) 
is the thermal velocity of X, and a gr is the grain cross section). The time-scale to is simply 
(tdeptcr) I '(tdep + t cr )- We are assuming here that one can neglect destruction of X due to either 
gas or solid-phase reactions. 

This rather crude model of depletion has the advantage that it can be rapidly computed as a 
function of density and time. For the present purpose, we (as in CWT99) continue the calculation 
until such time as we reach a central CO column density compatible with the observed C 17 
column density toward the central dust peak in L1544. At this point, as discussed by CWT99, 
the main contributions to the observed CO column come from lower density foreground and 
background material while CO in the central high density peak is depleted to very low abundances. 
The present calculations confirm this. Even though cosmic ray desorption is included, the central 
CO abundance is still depleted by a factor of 10 3 (~ 30), if the WMA (BAP) density structure 
is used. This also depends sensitively on Ed(CO) and indeed our results suggest that Ed(CO) 
must be larger than ~ 900 K in order to have depletion in the central region as observed (for the 



standard values of C and S, see Tab. H). A similar result has also been found by Aikawa et al. 



2001, who cannot reproduce the observed "hole structure" in the CO column density if the grain 



surface is covered by non-polar ice (i.e. -Ed (CO) = 960 K). 

Although the gas phase chemistry in general has long time scales, the ion chemistry is 
relatively rapid and thus one can expect in reasonable approximation that the abundances of ionic 
species such as HCO + , N2H+ etc. are given by the steady state chemical equations using the 
instantaneous abundances of the neutral species which are important heavy element repositories. 
The time scale for the "ion chemistry" is expected to be determined by recombination to species 
such as HCO + with a rate of Pdiss n ( e ) where Pdiss is the dissociative recombination rate and n(e) 
the electron density. For "canonical estimates" n(e) ~ 10~ 3 cm -3 and j3di ss ~ 10~ 6 cm 3 s _1 , we 
find a time scale of 30 years for the ion chemistry which is much less than depletion time scales. 
Thus one can compute HCO + in terms of the instantaneous CO abundance and N2H + in terms of 
the instantaneous N2 abundance. 

Analogously, the electron fraction x(e) can be computed in terms of global estimates for the 
molecular and metallic ions and using the same instantaneous abundances for CO, N2 etc. in 
the gas phase. The approach we have adopted is a simplified version of the reaction scheme of 
Umebayashi fe Nakano 199C| . Thus, we compute a generic abundance of molecular ions "mH + " 
assuming formation due to proton transfer with H3" and destruction by dissociative recombination 



with electrons and recombination on grain surfaces (using rates from Draine & Sutin 1987). The 



-12- 



abundance of Hg~ is calculated considering formation as a consequence of cosmic ray ionization of 
molecular hydrogen and destruction by proton transfer with heavy molecules (essentially CO and 
N2 here) as well as due to the processes mentioned above for molecular ions. We also consider 
metal ions M + formed by charge transfer and (including the equation for charge neutrality) solve 
a cubic equation for the electron abundance. 

In similar fashion, we have computed the instantaneous [DCO + ]/[HCO + ] (and [N2D + ]/[N2H + ]) 
ratios assuming that the respective deuterated species form by proton transfer from H2D + . Thus, 
we express the abundance ratio Rd = [DCO + ]/[HCO + ] as : 

R D = (1/3) ki x{HD) /%x(e) + k 2 x(m) + k 3 x(gr)). (5) 

In the above, k\ is the rate for production of H2D + due to reaction of HD with H3", k e is the rate 
for dissociative recombination of H2D + , hi is the rate coefficient for (and H^D" 1 ") destruction 
with neutral species, and k% is the rate for recombination of H2D + on grain surfaces ( praine 



fc Sutin 1987[ ). x(HD) is the HD abundance relative to H2 taken to be 3 x 10 5 ( [Linsky et al 



1995), x(m), a function of the amount of depletion, is a weighted average over neutral molecules 
which can undergo proton transfer with H^, and x(gr) is the grain abundance by number. The 
latter has been computed (as has k^) using a MRN QMathis et al. 1977 ) distribution with lower 



cut-off radius a m j n (100 A in standard case) and upper cut-off 0.25 fim. This treatment ignores 
effects due to atomic deuterium flDalgarno fc Lepp 1984ft and moreover assumes that all deuterium 
enrichment originates in H^D" 1 " (rather than e.g. CH2D + ). However, it improves on equation (1) 
of CWT98 in that it treats the depletion (via x{m)) in consistent fashion and in that it considers 
the effect of recombination on grains. From eqn. |5] it is clear that an increase of depletion boosts 
the deuterium fractionation. In fact, depletion causes a drop in the H2D + and H3" destruction 
rates, and a rise in the H2D + formation rate, via Hg~+ HD, due to the larger H3" abundance. 
The net result is a larger H2D + /H^ abundance ratio and a consequently more efficient deuterium 
fractionation. 

In Fig. ^, we show the dependence on radius of various ionic and molecular abundances 
predicted by our standard models fitting the central CO column density of L1544 (left panels 
assume the density structure fit to the results of WMA, whereas the right hand panels were 
obtained assuming the density structure deduced by BAP). The top panels show predicted 
depletion factors for CO and N2. The bottom panels shows the electron fraction, and the predicted 
abundances for the various ionic species observed by us. We stress that this is not our "best fit" 
model (see Sect. |3.2.2j) but it illustrates several features of our results. 

We note first the rapid increase of the CO depletion factor with decreasing radius, i.e. the 
gaseous CO abundance is reduced to a value of order 10~ 3 (0.03) of the "canonical value" of 10~ 4 
relative to H2 in the constant density central region, using WMA (BAP) density values. This 
contrasts with N2 which (although itself depleted) becomes the most abundant heavy molecule 
in the center. It is the rapid variation of these quantities as a function of density which renders 
necessary a consideration of abundance variations along the line of sight and hence radially. Also 
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shown in Fig. ||] are the variations of the DCO + and HCO + abundances. Their ratio (Rd) increases 
by roughly an order of magnitude from 0.04 (0.04) to 0.7 (0.35), using WMA (BAP) parameters, 
between the edge and center of our model. 

The variation in ion densities shown in the bottom panel of Fig. |I| reflects the behavior of the 
depleted species in the two standard models. Thus, in a WMA core, HCO + abundance is strongly 
reduced in the center because of CO depletion, whereas DCO + shows a rather flat dependence 
due to the increase of the deuterium fractionation with increasing /r>. A similar behaviour is 
present for the N2H + and N2D + abundances, with the consequence that Rd and [N2D + ]/[N2H + ] 
become greater than 1 in the innermost parts of the core. In a BAP core, the lower density 
causes less efficient CO and N2 depletion, and much shallower profiles for HCO + and N2H + 
fractional abundances. The ions plotted in the figure are "major ions" at all radii (except that 
becomes abundant in the central region under some circumstances) suggesting that for depleted 
high density cores, the sum of the observed molecular ion abundances may give a good estimate 
of the electron abundance. This is a great simplification relative to "normal cores" where "metal 
ions" have comparatively high abundances. Our calculations suggest that in high density depleted 
cores, the depletion of "metals" (assumed here to behave like CO) is likely to be so large that 
this complication becomes negligible. We have used an initial metal abundance about one order 
of magnitude lower than that measured in diffuse clouds, as is usual ("low metals") in chemical 



models of dense clouds (e.g. |Prasad fc Huntress 1982[ |Herbst fc Leung 1989| ; |Graedel et al. 1982 



Lee et al. 1996|) . 



The most convincing aspect of our model is that it qualitatively explains the enhancement 
of N2D + in the center of L1544, close to the depletion peak. Thus it explains the higher degree 
of fractionation of deuterium in N2D + as compared with HCO + . It also explains roughly the 
differing depletions observed for N2H + , CO, and HCO + . In fact, as discussed in the next section, 
the qualitative features of the model are in excellent agreement with observation. The pattern 
that ions and deuterated species peak close to dust emission maxima appears to be a general 
feature (Tafalla et al., in prep.) of nearby cores and suggests that the model presented here has 



wider applications than to L1544. 



3.2.2. The best fit model 

Several model calculations have been performed to find the "best fit model", i.e. the model 
which best reproduces observed column densities at the peak of L1544. We started with a 
"standard" cloud model, where Q = 1.3xl0 -17 s _1 , the CO binding energy E'd(CO) = 1210 K 
( Hasegawa et al. 1992 ), the N2 binding energy E-£>{^2) = 787 K, the minimum radius of dust 

grains ttrnin 

= 10 6 cm, the sticking coefficient S = 1.0 (see Table ||). Two density distributions 
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(from WMA and BAP) were considered and the above parameters were varied until the quantity: 

-i 2 



x 2 = E 



=1 L 



' (JVobs(») - N mod (i)) 



(6) 

a N obs (i) 

was minimised. In eqn.(^), N ^ s (i) and N mo( i(i) are the observed and model calculated column 
density of HCO+, DCO+, N 2 H+, and N 2 D+ at the L1544 "molecular peak". a Nohs{i ) is the 
uncertainty associated with N Q i, s (i), which has been assumed here to be a 30% error. Thus we 
assume the errors in our column density determination to be dominated by systematic effects such 
as those caused by our summary treatment of the radiation transfer. 

The column density dependence as a function of impact parameter b predicted by the two 
"best fit" models (for the two density distributions) can then be compared with those observed. 
We have made azimuthal averages of our observed column densities using averages within bins 
defined by i < b < i + 20 arcsec, with i = 0, 20, 40, ... (arcsec), with the exception of the value at b 
= 0. This is shown in Fig. |5[ where large symbols represent the average column density inside the 
corresponding bin (symbols mark the upper edge of each bin). One sees that there is considerable 
scatter as indeed one can expect given the elongated nature of the maps. This is in particular the 
case for N2H + where the "mean" calculated in this fashion may clearly be misleading. One sees 
nevertheless that the distributions become clearly more peaked going from CO (flat) to N2D + . 

In Fig. ||, the observed column density profiles are compared with those predicted by four 
models compatible with the peak column densities. The parameters used in these models are 
reported in Tab. ^together with the x~squared estimates using equation @. To understand our 
results, it is useful to recall that we need to assume high degrees of CO depletion in order to 
explain the CO and HCO + column densities. A certain amount of nitrogen depletion is also needed 
to explain N(N2H + ). Such low gas phase molecular abundances cause extremely high degrees of D 
fractionation however with the model expectation for e.g. [DCO + ]/[HCO + ] higher than unity in 
the dense depleted central regions. This then gives difficulty explaining the abundances of N 2 D + 
and N2H + . The range of parameters in Tab. [2| demonstrates that the peak column densities on 
their own only give partial constraints on the chemical models. 

Fig. ^ however shows that the extent of the emission in the various ions is also an important 
constraint. In models 1 and 2 (Tab. |2j) for example, the model half-width for N2H + is much too 
small. While this could be partially due to our azimuthal averaging (see Fig. ||), it is also the case 
that the cosmic ray ionization rate has been chosen to have values in models 1 and 2 more than 
an order of magnitude lower than in the standard model and the result is that N(N2H + ) is much 
too small at large radii. 

A way of improving the agreement between observations and model predictions is to more 
efficiently deplete CO outside the core center (to increase the N2H + column density at b > 0). 
This can be done by increasing the CO binding energy and thus reducing the efficiency of cosmic 
ray induced desorption. However, this also causes an increase in deuterium fractionation which is 
not observed. As discussed in Sect. 2.4.1| , [DCO + ]/[HCO + ] is observed to have rather low values 
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at high impact parameter. This suggests that reactions other than with CO and N2 act to reduce 
the H2D + abundance. We therefore postulate the existence of a volatile neutral species which will 
remain in the gas phase when CO is depleted and which has a proton affinity permitting transfer 
of a proton from H^. 

This species needs to be abundant (at least as abundant as CO) to effectively keep the 
deuterium fractionation low in those regions where CO is significantly depleted. One possibility 



may be atomic oxygen, which is predicted to be quite abundant in dense clouds (e.g. [Lee et al 



1996 ; Bcrgin et al. 2000 ; Viti et al. 200 1|) . Its binding energy onto a dust grain covered with 



water ice is thought to be a factor of about 1.5 smaller than that of CO (-Ed(O) ~ 800 K; Tielens 
&; Allamandola 1987). We have therefore run several models which include atomic oxygen in the 
chemistry with different values of -Eq(O), and determine the x 2 an d column density profiles in the 
same fashion as before. These results are shown in Fig. |7| and correspond to models 3 and 4 of 
Table ^. One sees that model 3 in particular gives a reasonable fit although the DCO + column 
density profile is still somewhat deviant. Both models 3 and 4 require an oxygen binding energy 
of ~ 600 K. 

Using the "best fit" models in Tab. namely Models 1 and 3, we determined the variation 
of the electron fraction with cloud radius r. Fig. |8| shows the fractional abundance of electrons, 
HCO+, DCO+, N 2 H+, and N 2 D+, as well as the depletion factor / D , as a function of r. Model 3 
(the best fit model) predicts a [N2D + ]/[N2H + ] aboundance ratio equal to ~ 0.4 at radius r ~ 2500 
AU. Both models show x(e) between ~ 10~ 9 at r ~ 2500 AU (n(H 2 ) ~ 10 6 cm" 3 ) and ~ 10" 8 
at r ~ 10 4 AU (n(H2) = a few x 10 4 cm -3 ). In these models, the dependence of x(e) upon gas 
density n(H 2 ) is approximately given by : 



x(e) = 6.3 x 10~ 6 x n(H 2 ) 

in Model 1, and 

x(e) = 5.2 x 10~ 6 x n(H 2 ) 



-0.64 



-0.56 



in Model 3. These estimates are roughly an order of magnitude lower than the standard relation 
between x(e) and ra(H 2 ) (x(e) = 1.3xl0 -5 x n(H_2)~ ' 5 ; McKee 1989| ). In other words, our results 



are compatible with electron fractions only a factor of order 3 higher than the lower limits based 
on the inferred HCO + abundance in the nucleus of L1544. Such low electron fractions are linked 
to the depletion of metals which, being positively charged, directly affect the electron budget of 
cloud cores. We use an initial abundance of metals 10 times lower than in the standard relation 



quoted above and this causes the difference in the coefficient (see Sect. 3.2.1 ). The difference in the 
exponent is to be attributed to the further metal depletion inside the core. If larger initial metal 
abundances are used, x(e) increases at the outer edges of the core, but it does not significantly 
change at the core center, where metal abundances are negligible because of depletion. In this 
case, a steeper slope in the x(e) — n(H2) relation is obtained. We caution that these estimates 
have large uncertainty as the differences in the models of Tab. [2] shows. Our approach to inferring 
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molecular ion densities as a function of radius is questionable and the real geometry is non 
spherically symmetric. There are additionally many uncertainties in the chemistry. Nevertheless, 
we conclude tentatively that the "standard relationship" may give a considerable over-estimate of 
the electron fraction in L1544. This has consequences for the ambipolar diffusion timescale. 

In the central portion of L1544, ionization degrees of the order of 10 -9 imply an ambipolar 
diffusion time scale of roughly 2.5xl0 13 x{e) = 2.5x10 yrs (see Spitzer 1978 and |Shu et al. 1987 ), 
comparable to the free-fall time scale at densities of ~ 10 6 cm -3 . This is consistent with the view 
that within the central region, (see |Ciolek and Basu 2001 , Ciolek and Basu 2000| ), conditions 
are "super-critical" and the core is rapidly developing towards a situation where it will collapse. 
However, the "caveats" noted above mean that one needs to test this conclusion with observations 
capable of better delineating the structure of the high density nucleus of L1544. 



3.2.3. Chemical problems with the best fit model 

There are some problems related with the existence of atomic oxygen in the gas phase which 
need to be clarified. In our models, the large abundance of gaseous O is due to the low value of 
£"d(0) and to the cosmic ray desorption mechanism which allows a prompt return of this species 
to the gas phase, once adsorbed onto grain surfaces. We are neglecting here the fact that the 
reactions of atomic oxygen with Hjj" necessary to reduce deuterium fractionation will also produce 
water and OH in the gas phase which can then deplete out onto grains, although only a small 



fraction of O (< 0.01) will be lost due to this process (e.g. Lee et al. 1996). An important 
consequence of the Hg~+ O reaction is that H30 + becomes the most abundant molecular ion in 



the depleted core (x(H 3 0+) ~ 10~ 9 at ra(H 2 ) ~ 10 b cm as also found by |Aikawa et al. 200l[) . 
However, the destruction rate of atomic oxygen through reactions with H3" ( r n+) is significantly 
slower than the O accretion rate onto dust grains (r^ust) []• This means that if the accreted O will 
not be completely processed onto grain surfaces, through, say, hydrogenation leading to water, 
between successive cosmic ray bombardments, we expect a fraction of free oxygen to be maintained 
on the gas phase. Proving the validity of these assumptions is outside the scope of this paper but 
we consider briefly the likely chemistry of atomic O on grain surfaces. 



Following Hasegawa fc Herbst 199^ ( see [Leger et al. 1985[ ), if the major source of nonthermal 



grain heating is due to Fe nuclei with energies 20-70 MeV nucleon -1 , the time interval between 
successive cosmic ray impacts is about 10 6 years. If the fractional abundance of atomic oxygen 
is ~ 10~ 4 and n(H) ~ 1 cm 3 , as gas phase models of dense clouds predict at steady state for 



2 At a density of ~ 10 6 - 10 5 cm -3 , our best fit model predicts x(Ut) ~ 1 - 2 x 10 10 , so that r„+ = 

3 

k„+ x xCRf) x n(H2) ~ 1 - 2 x 10~ 14 cm _3 s _1 , where the rate coefficient k„+ ~ 10~ 9 cm 3 s -1 . For the O 

3 3 

accretion rate, assuming a unity sticking coefficient and dust grains of 10 -5 cm in size, r-d us t ~ 10~ 17 n(H.2) cm -3 , 
implying rdust/r H + ~ 50 - 100 at n(H2) = 10 5 - 10 6 cm -3 , respectively. 
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£ = 1CT 17 s- 1 (e.g. |Lee et al. 19961 ), T = 10 K, and n(H 2 ) ~ 10 5 cm" 3 , the O and H accretion 
rates are ~ 4xl0~ 5 s _1 and lxl0~ 5 s -1 , respectively, i.e. between two cosmic ray bombardments 
about 10 9 O and 3 x10 s H atoms can accrete on the surface of a grain. Assuming that hydrogen 



can quickly move on the grain surface (from the laboratory work of Katz et al. 1999, the H 
diffusion rate on silicate grains is about 2xl0 -4 s , much larger than the c.r. heating rate), a 
fraction between 15% and 30% of surface O will be hydrogenated and transformed in OH and 
H2O. The rest of the O atoms may form O2, but observations have shown that this is probably 
not the main repository of oxygen either in the solid QVandenbussche et al. 1999 ) or in the gas 
phase ( poldsmith et al. 2000 ). It is thus possible that a large fraction of atomic oxygen remains 
unreacted on the grain before the next cosmic ray will release it back in the gas phase. 

It is relevant in this contest that the detailed gas-grain chemical-dynamical model of L1544, 
presented by Aikawa et al. 2001| (which includes depletion onto dust grains but no surface 
processing of accreted species) predicts large O fractions across the core (from ~ 10~ 4 at n(H2) 
~ 10 cm " 3 to a few x 10~ 5 at 10 6 cm , similar to our results). This is in their "fast collapse" 
model, which fits the observed CCS and CO distributions. Our model is crude compared with 



that of Aikawa et al. 2001, but nevertheless both models reproduce the observed features quite 



well, proving the validity of our chemical assumptions in Sect. |3.2.1 . There are some differences 



between the two models concerning the prediction of N2H + and N2D + abundances. In Aikawa et 



al. 2001| , the calculated column densities of the two species are smaller than those observed by 
about an order of magnitude, perhaps due to uncertainties in the N2 formation mechanism. We 
"avoid" this problem by simply starting the chemistry with all the available nitrogen already in 
the form of gas phase N2 (Sect. 3.2. 1| ), thus reproducing the observed N2H + and N2D + column 
density variations (see Model 3 in Fig. 0). 



Conclusions 



An important result of this study is that N2H + shows more deuterium fractionation 
towards the dust emission peak of L1544 than HCO + (models without depletion gradients expect 
[DCO+]/[HCO+] = [N 2 D + ]/[N 2 H+] [Rogers and Charnley 2001 ). We suspect that the high degree 
of deuteration observed in ammonia (see [Tine et al. 2000] , |Roueff et al. 2000|) m some cores can 
be best explained by a model similar to that which we have adopted. Ammonia like N2H + is 
easily produced from N2 (see Rogers and Charnley 200l| ) and is therefore likely to be relatively 
abundant in the depleted region. The compact nature of "ammonia cores" is naturally explained 



if this is the case ( Tafalla et al., in prep. ) 



We have attempted to obtain density estimates for the regions within the L1544 core giving 
rise to the observed DCO + and N2H + (N2D + ) emission. We obtain in this way values between 10 5 
and 10 6 cm -3 consistent with the general idea that we are observing selectively the high density 
layers. The results from DCO + suggest densities as high as 10 6 cm -3 consistent with the values 
inferred from dust emission and absorption. However, lower values are obtained based upon our 
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N2H + and N2D + data. This is surprising and we are presently unsure of the explanation. The 
observed spatial distribution suggests that N2H + and N2D + reside in higher density layers than 
DCO + . The reasons for this discrepancy could have to do with the isothermal assumption which 
we have made based upon our C 18 and NH3 results. Attempts to simulate the dust emission from 
L1544 ( ^ucconi et al. 2001 , Evans et al. 2001] ) have shown that there is probably a subtantial 



fall off in dust temperature from the center to the edge with values of order 7 K in the nucleus. 
At densities above 10 5 cm~ 3 ( [Krugel and Walmsley 1984 ) gas and dust temperatures are probably 



coupled and hence there may be a gradient in the gas temperature too (see also Goldsmith 2001 ). 
This should be taken into account in future studies. 

We have developed a crude model of the ion-chemistry in the core of L1544. This simulates 
the observed depletion and can reproduce the observed dependence of the column densities of 
species such as N2H + ,N2D + , HCO + , and DCO + as a function of offset. The main advantage of this 
is that it allows us in objective fashion to consider the relative contributions to observed column 
densities from the high density depleted nucleus and the lower density foregrond (background) gas. 
There are large uncertainties in both the input to the chemistry and the process of inferring radial 
dependences of molecular abundances from the observations. One interesting result is that we get 
the best fit to the observed deuterium fractionation in models in which atomic oxygen is allowed 
to remain with an abundance of order 10 -4 in the gas phase. There are some problems associated 
with this but one needs more detailed models of both the surface and gas phase chemistry to test 
this. It does have an observational consequence which is that H^O" 1 " becomes a major ion towards 



the dust peak, in agreement with numerical chemical models of Aikawa et al. 2001 



Applying this model, we find that the ionization degree in L1544 is likely to be an order of 
magnitude smaller than estimated using "canonical formulae" existing in the literature, mainly 
because of the reduced metal abundances. The result is a tentative finding because we have 
difficulty in fitting the observed abundances as a function of offset in LI 544 and in particular 
the observed degree of D fractionation. However, our estimates for ionization degree do not vary 
greatly when one compares the different models which give an adequate fit to the observed column 
density distributions. We conclude that in the case of L1544, the ionization degree is lower and 
hence the ambipolar diffusion timescale is shorter. Our estimated timescale is of the same order 
as the free-fall time consistent with the idea that the nucleus of the LI 544 core is undergoing 
dynamical collapse. 

A final point to emphasize is that in models without atomic oxygen, or where O also 
experiences strong depletion, the main ion in the highly depleted region is H^, so that the 
mean mass of positive ions decreases from ~ 30 (in the "canonical" undepleted case, where the 
representative ion is HCO + ) to ~ 3. This has the consequence of further reducing the ambipolar 
time scale, due to a drop in the drag coefficient, by a factor of 7-8 (D. Galli, priv. comm.). 

Other cores may yield different results and the effects of the possible temperature gradient 
mentioned above need to be taken into account. Finally, it will be important in future studies to 
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take into account a more realistic geometry (non spherically symmetric) than we have done. 
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A. Column density determination 

In addition to the LVG calculations, we also used the following expressions to determine the 
total column density. For optically thick transitions: 

Ar TOT = 8 ^ /2Av 9L T _ 9ML , (Al) 

2VTn2X 3 Ag u 1 - exp(-hu/kT cx ) giexp(-Ei/kT ex ) 

where Av is the line width, A and v are the wavelength and frequency of the observed transition, 
respectively, A is the Einstein coefficient, g\ and g u are the statistical weight of the lower and upper 
levels, r is the optical depth, h is Planck constant, T ex is the excitation temperature (assumed the 
same for all rotational levels), Qrot is the partition function, E\ is the energy of the lower level, 
and k is Boltzmann constant. For linear molecules (as those observed in this paper): 

oo 

Qrot = E( 2J + l)exp[-Ej/(KT)] (A2) 
j=o 

Ej = J(J+l)hB, (A3) 

where J is the rotational quantum number, and B is the rotational constant. For rotational 
transitions with hyperfine structure (e.g. N2H + (l-0) and N2D + (2-l)), r refers to the total optical 
depth (given by the sum of the peak optical depths of all the hyperfine components), and Av to 
the intrinsic line width. The error on Atot is given by propagating the errors on Av, r, and T ex 



in eqn. Al 



If a line is optically thin and all the rotational levels are characterized by the same excitation 
temperature T ex , the expression of the total column density becomes: 



iVTOT = - - ^ , (A4) 

X 3 A g u J u (T ex ) - Ju(T hg ) 1 - exp[-hu/(kT ex )] giexp[-Ei/(kT eyL )] 

where W is the integrated intensity of the line (W = \JH j {2\J InT) x AvT m b, for a gaussian 
line, with T m b = main beam brightness temperature), J u (T ex ) and J u (T\> g ) are the equivalent 
Rayleigh-Jeans excitation and background temperatures. The integrated intensity of each line has 
been measured by integrating over a velocity range determined in the following way: (i) sum all 
the spectra in the map; (ii) find the rms, or the 1 a level of the noise in the off-line channels of 
the sum spectrum; (iii) include in the velocity range for integration (at zero level) all the channels 
in the line which are above the 1 a level determined in point (ii). For the calculation of the total 
column density only those positions with W/aw > 3 have been considered {aw is the error on W, 
determined by the expression aw = Av res x rms x \/N c h, where Au res is the spectral resolution 
in kms" 1 , and iV c h is the number of channels in the integrated area). The error on A^tot is simply 
given by a NTOT = a w x N TOT /W, 

The parameters used to determine Atot for each species are listed in Table |l[ 
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A.l. DCO+ and HCO+ 



As we saw in Sect. 2.2, an LVG statistical equilibrium code has been used to estimate number 
density and DCO + column density from observed DCO + (3-2) and (2-1) lines. To quantify the 
uncertainties in column density estimates from statistical equilibrium (SE) calculations, due to 
the unknown DCO + fractional abundance and H2 number density, we first determined DCO + 
column densities from the SE code and a least square method applied to model calculated data 
points. To do this, the SE program was run for different values of H2 number density (between 10 3 
and 10 6 cm -3 ) and DCO + fractional abundances (between 10 -12 and 10 -9 ). Each run furnishes 
a value of T ex , r, T m t, of the first four rotational transitions of DCO + , together with the DCO + 
column density divided by line width. Secondly, a least-square method was applied to these model 
data points to find general expressions for estimating iV(DCO + ) once the brightness temperature 
T mh{2 ~i) of the J = 2-1 line and the iy(DCO + (3 - 2))/W(DCO+(2 - 1)) (= R w ) integrated 
intensity ratio are known from observations. However, this gives -/V(DCO + ) values with associated 
errors greater than 30% (N/a^ < 3), so we decided to first determine the excitation temperature 
(T ex ( 2 -i)) and the optical depth (t( 2 -i)) of the J = 2-1 line using the least-square method, which 
gives: 

r ox(2-i) = 01 + a 2 logR w + a 3 T mb (2 - 1) (A5) 
logr^ = b 1 + b 2 (logR w ) 2 + b 3 log{T mh (2 - 1)), (A6) 

with 01 = 11.3±0.1, a 2 =10.3±0.2, a 3 =-0.05±0.01, &i=-0.2±0.1, 6 2 =2.0±0.5, &3=1.0±0.1. 
Secondly, we assumed T ex = T ex ( 2 _i) for all rotational levels, and finally calculate N using 
expressions (|A3J) or ( |Al| ) in the appendix, depending on the value of the optical depth (t( 2 _i) 
< or > 0.5, respectively). Only positions inside the half maximum contour of the DCO + (2-l) 
integrated intensity map have been included in this calculation, to exclude low sensitivity spectra 
from the analysis. The particular form of these expressions was chosen because of the relatively 
low x 2 - If from ( A6) t/o> < 3, TY2-1) was estimated from the radiative transfer equation 



r( 2 -i) 

-^mb 



-log 



1 



where J u (T ex ) and J v (Tbg) are the equivalent Rayleigh- Jeans excitation and background 



(A7) 



temperatures, and T ex is given by eqn. (|A5|) . Finally, if the resultant A^/ctn < 3, T ex was fixed to 7 
K, its mean valued], and r estimated from ( |A7[ ). T cx = 7 K was also assumed for all the positions 
outside the half maximum contour of the DCO + (2-l) integrated intensity map, where eqns. (|A5|) 
and ( |A6| ) cannot be applied because of the poor i?vv estimate. In the particular case of the (20, 
-20) offset position, the optically thin D 13 CO (2-1) line was used to estimate the DCO + column 
density, assuming the same excitation temperature found for the DCO + (2-l) line. We note that 



3 T ex drops from ~ 9 K towards the map peak, at offset (20, -20), and two adjacent positions to < 7 K at larger 
distances from the peak. 
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the different approaches used in estimating column densities do not change the results by more than 
a factor of 2. However, it is important to estimate the optical depth of the line before determining 
the column density. For example, in the case of offset (20, -20) where the DCO + (2-l) line has one 
of the largest values of optical depth (~ 2), (i) the SE calculation together with equation (Al) 



give iV(DCO + ) = (4±l)xl0 12 cm" 2 , (ii) the use of T cx = 7 K and eqn.(|57|), gives r = 1.7±0.2 
and iV(DCO+) = (2.0±0.2)xl0 12 cm" 2 , and (iii) the use of D 13 CO + (2-l) leads to iV(DCO+) = 
(4.0±0.3) x 10 12 cm -2 . If the optical depth of the DCO + (2-l) line is not taken into account, one 
obtains iV(DCO + ) = (9.6±0.2)-(10.2±0.2) x 10 11 cm" 2 , if T ex = 7 - 9 K, respectively. 

For the determination of the HCO + column density, the excitation temperature was assumed 
equal to that found for DCO + . We used the optically thin HC 18 O + (l-0) line, if available, and 
the H 13 CO + (1-0) line in all other positions. The passage from 7V(HC 18 + ) and iV(H 13 CO + ) to 
iV(HCO + ) was made by multiplying the former quantities by 560 and 77, respectively, based on 
the [ 16 0]/[ 18 0] and [ 12 C]/[ 13 C] abundance ratios in the local ISM flWilson fe Rood 1994j) . 

The optical depth of H 13 CO + (1-0) was calculated from eqn. (|A7|). As in the case of DCO + , 



eqn.(A4) or (Al) was used, depending on the value of r. It is interesting to note that the resultant 
r in many positions is less than 0.5, i.e. H 13 CO (1-0) is mostly optically thin and thinner 
than the DCO + (2-l) line in the central part of the core. However, the profile of H 13 CO + (1-0) 
suggests the existence of absorption by a low density foreground layer (see Tafalla et al. 1998] , 



Williams et al. 1999[) which may cause a decrease in the brightness temperature and, consequently, 



in the line optical depth (see eqn. |A7|) . The three hyperfine components of the HC 17 (1-0) 



transition detected at offset (20, -20), and shown in Pore et al. 2001 together with their recent 
measurement in the laboratory, allowed us to check the thickness of the HC 18 O + (l-0) line. 
Assuming T ex = 9 K, the excitation temperature found from DCO + data at the same position, 
we find Ar(HC 18 + )/N(HC 17 + ) = 4.2±0.4, quite close to the [ 18 0]/[ 17 0] ratio (= 3.65; Penzias 



1981) expected in optically thin conditions. 



The determination of the HCO + column density starting from H 13 CO is affected by an 
addictional complication, i.e. the presence of 13 C fractionation in cold and dense gas (e.g. Smith 



fe Adams 1984) due to the exothermic reaction: 

13 C++ 12 C0 _^ 12 C + + 13 CO + AE, 

with AE/k = 35 K (Watson 1977). In optically thin conditions, the A^(H 13 CO + )/A^(HC 18 + ) 
column density ratio should be equal to the product of [ 13 C]/[ 12 C] and [ 16 0]/[ 18 0] abundance 
ratios (~ 7). In L1544, N(H 13 CO + )/N(HC 18 + ) ranges from 3 (at offset [80, -80]) to 7 (at 
offset [40, -80]) with an average value of 4dbl. It is thus always lower than the local interstellar 
medium value, suggesting that 13 C fractionation is not significantly affecting our conclusions. 
However, given that H 13 CO + (1-0) is probably affected by foreground absorption, it is extremely 
hard to estimate the effects of 13 C fractionation. From the current data we determined the 
H / [H 13 CO + (l — 0)]/VF[HC 18 O + (l — 0)] integrated intensity ratio and found an average value of 
4±1, smaller than the value expected in optically thin conditions and no 13 C fractionation. This 
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result suggests that optical depth effects on the H 13 CO (1-0) line are in any case predominant. 



A.2. CO 

For the CO column density, we used C 17 O(l-0), when available, and C 18 O(l-0) in the other 
positions assuming T cx = 10 K (Sect. [Hp, [ 17 0]/[ 16 0] = 2044, and [ 17 0]/[ 16 0] = 560 ( [Wilson 



& Rood 1994, Penzias 1981| ). C 17 O(l-0) has (well resolved) hyperfme components with relative 
intensity ratios consistent with optically thin emission, whereas C 18 O(l-0) has a moderate optical 
depth (r between 0.5 and 1), based on the comparison between the W / (C 18 0(1 — 0))/W(C 17 O(l — 0)) 
integrated intensity ratio and the [ 18 0]/[ 17 0] abundance ratio in the local ISM. Following [Myer^ 



et al. 1983| , we derive: 

T mb [C 18 0(l - 0)] _ „ „ A - exp(-T 18 ] 



3.65 x ^ — i5i (A8) 

T mb [C 17 O(l-0)] \ r w J' K J 

where T m b(i) is the main beam brightness temperature of line i and T%g is the C 18 O(l-0) optical 
depth. Once ns is known from the above equation, the excitation temperature of the C 18 line 
can be determined from the radiative transfer equation: 

hv/k 



(A9) 



where v is the frequency of the C 18 O(l-0) line. The variation in from point to point in the 
map is not significant, considering the errors, and we assumed a constant value of T ex (= 10 K) 
given by the weighted mean of the excitation temperature in the observed positions. As in the 
case of A^(DCO + ) and A r (H 13 CO + ), when r > 0.5, the C 18 column density has been calculated 
using eqn. ( |Al[ ). 



A.3. N 2 H+ and N 2 D+ 

The observed N2H + (l-0) and N2D + (2-l) lines present hyperfme structure due to the 
interaction between the molecular electric field gradient and the electric quadrupole moments 
of the two nitrogen nuclei. The J = 1— »0 line of N2H + is splitted in seven components (see 
Caselli et al. 199~5[ ), whereas N 2 D + (2-l) has 38 hyperfines which partially overlap because of line 
broadening. 

We first determined the intrinsic linewidth, total optical depth, and excitation temperature 
from the hfs fitting procedure applied to high sensitivity spectra (selected by requiring W /ay/ > 
20, W being the intensity integrated under the seven hyperfme components). In fact, the optical 
depth determination becomes more and more uncertain with increasing spectral rms, which can 
alter the relative intensities of the components and thus crucially affect r estimates. Therefore it 
is important to have clear detections of the seven components before attempting an hfs fit. 
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However, even in the presence of high sensitivity spectra one can encounter problems if the 
transition is very thick (r > 30) because the CLASS fit procedure is limited to r values smaller 
than 30. In these cases, it is convenient to estimate the N2H + (l-0) optical depth from the thinnest 
component (the lowest frequency Jfi,f = li,o ~~ ► 0l,i nne > see paselli et al. 1995 ) assuming a 



certain value of T ex (e.g. the mean value, < T ex >, which in L1544 is 5 K), and using eqn. A.7. 
The hfs fit gives r = 30 - with a suspiciously low associated error (~ 0.1) - in four positions of 
the L1544 N 2 H+(l-0) map (offsets [40, -40], [20, -40], [40, -60], and [20, -60]). For these spectra, 
the use of eqn. with T cx = 5 K confirms the presence of high optical depth (< 30). In these 
cases, we used the integrated intensity of the "weak" and moderately thick (7-1,0-1,1 < 1) = 
1,0 — ► 1,1 component to determine the total N2H + column density from eqn.( |A4[ ). Although this 
may underestimate the total column density by a factor of ~ r/(l — exp(—rl, — 1, 1)) < 1.6, 
the uncertainty associated on the above estimate of the total optical depth is also large (about 
20%), and does not take into account of the possible presence of excitation anomalies, as found by 



Caselli et al. 1995[ For consistency, the use of the Fi,F = 1,0 — > 1,1 component to estimate the 



total N2H + column density has been extended to all other positions where this hyperfine is well 
detected. 

For those N2H + (l-0) spectra which cannot be hfs fitted because of low sensitivity, and where 
the "weak" component is not visible, we assumed optically thin conditions, T ex = 5 K, and 



determined A r (N2H + ) using eqn. A4 



As we already pointed out at the beginning of this section, the J = 2 — ► 1 transition of 
N2D + is splitted in 38 hyperfine components and only in one position (the integrated intensity map 
peak, where the N2D + (2-l) observations have been repeated several times to check the system) 
it has been possible to determine the total optical depth and excitation temperature from the 
hfs fit (r to t = 4.9±0.6, T ex = 4.9±0.8 K). In all the other positions we (i) assumed optically thin 
conditions, (ii) estimated the integrated intensity below the hyperfines (W) with its uncertainty 
(<7w)j (ih) excluded all spectra with W /aw < 3, (iv) fixed T cx at 5 K, and (v) use eqn. (A4) to 
estimate the total column density. 
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Table 1: Molecular parameters used to calculate iVroT 



Transition 


B 




Av rcs 


iVch 




MHz 


Debye 


krns" 1 




N 2 D+(2-l) 


38554.719 


3.4 


0.038 


275 


HC 18 O + (b-0) 


42581.21 


3.9 


0.034 


17 


H 13 CO + (1-0) 


43377.32 


3.9 


0.034 


25 


HC 17 O + (l-0) 


43528.933 


3.9 


0.067 


6 


N 2 H+(l-0) 


46586.867 


3.4 


0.063 


75 


C 18 O(l-0) 


54891.420 


0.11 


0.027 


33 


C 17 O(l-0) 


56179.990 


0.11 


0.026 


72 


D 13 CO + (2-l) 


35366.712 


3.9 


0.041 


14 


DCO + (2-l) 


36019.76 


3.9 


0.041 


22 


DCO+(3-2) 


36019.76 


3.9 


0.054 


13 
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Table 2. Parameters of the "best fit" models 



Parameters Standard Model 1 Model 2 Model 3 Model 4 



Density Dist. 


WMA 


WMA 


BAP 


WMA 


BAP 


C(s- X ) 


1.3xl0~ 17 


5.0xl0~ 19 


1.3xl0~ 18 


6.0xl0~ 18 


5.0xl0~ 18 


£ D (CO)[K] 


1210 


960 


1210 


1210 


1210 


£d(N 2 )[K] 


787 


600 


600 


787 


787 


E D (0)[K] 








600 


600 


a m in(cm) 


lxl(T 6 


2.5xl0~ 5 


l.OxlO" 5 


5.0xl0~ 6 


1.0xl0~ 6 


S 


1.0 


0.1 


0.1 


1.0 


0.1 


x 2 




4.7 


6.7 


1.6 


5.4 
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Fig. 1. — Logarithmic values of the H2 number density, n(H.2), superposed upon the half maximum 
contours of the DCO + (2-l) (thin curve), the DCO + (3-2) (dashed curve) integrated intensity maps, 
and the 1.3mm continuum dust emission (dotted curve) convolved to a 22 " beam (circle in the 
bottom right). The grey scale shows the DCO + (3-2) integrated intensity map (levels 50, 70, and 
90% of the peak, 0.84 K kms _1 at offset [20, -20]). The density values in the figures are based on 
statistical equilibrium calculations (see text). 

Fig. 2. — Map of the depletion factor /d (grey scale) superposed to the 1.3mm continuum dust 
emission map of WMA, smoothed at a resolution of 22" (dotted contours), /d ranges between 1 
(at the core edges) and 10.8 (at offset [24, -34]); the grey scale (from light to dark) shows parts in 
the core with /d = 1, 3, 5, 7, and 9. Dotted contours levels represent 30, 50, 70, and 90% of the 
1.3mm map peak (225 mJy/22"at offset [26, -21]). The distance between the two peaks is about 
half a beam. 

Fig. 3. — (top left) DCO + column density map; contour levels are between 30 and 90%, in steps 
of 20% of the map peak (= 4.0 xlO 12 cm -2 , at [20, -20]). (center left) HCO + column density map; 
contour levels are as in the previous panel and the map peak is at offset (0, 0) (1.1x10 cm -2 ). 
(bottom left) iV(DCO + )/iV(HCO + ) column density ratio map. Contour levels are 10, 30, and 50% 
of the map peak (0.06±0.02 at offsets [20, 20] and [40, -20]). The difference between the four 
positions inside the 50% contour level is not significant, taking into account the associated errors. 
(top right) N2D + column density map (contours range between 30 and 90%, in steps of 20%, of the 
map peak (= 4.4xl0 12 cm -2 , at [20, -20]). (center right) N2H + column density map; contour levels 
are as in the previous panel and the map peak is at offset (20, -10) (2.0xl0 13 cm" 2 ), (bottom right) 
N(N 2 B + )/N(N 2 B. + ) column density ratio map; contour levels are 10, 30, 50, 70% of the peak (0.26 
at [40, -20]). 

Fig. 4. — Dependence on radius of the depletion factor for CO (/d) and N2 (top panels), electron 
fraction and various ionic species (bottom panels) for the "standard model" (see Tab. [2] for input 
parameters) using the density distribution found by WMA (left panels) and BAP (right panels). 

Fig. 5. — Column density (filled circles) versus impact parameter for (from top to bottom): 
CO, HCO + , DCO + , N2H + , and N2D + . White symbols are average column densities inside the 
corresponding bin, with the only exception of the value at b = (see text). 

Fig. 6. — Column density profiles predicted by the "best fit" models (solid curves) and deduced 
from observations (dotted curves). Left panels refer to models with a volume density profile from 
WMA (Model 1), whereas models shown in right panels assume the L1544 density profile from 
BAP (Model 2). 

Fig. 7. — Column density profiles predicted by the "best fit" models with atomic oxygen (solid 
curves) and observed (dotted curves). As in Fig. ||, left panels refer to models with volume density 
profile from WMA (Model 3), whereas right panels show models with density profiles from BAP 
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(Model 4). The inclusion of atomic oxygen in the gas phase furnishes a better agreement with 
observed molecular ion column densities at the dust peak and with observed column density profiles. 

Fig. 8. — Radial profiles of fractional abundances of electrons, HCO + , DCO + , N2H" 1 ", and N2D + 
predicted by Model 1 (top) and 3 (bottom; see Tab. |2[). Dashed curves represent the depletion 
factor. Electron fractions are about one order of magnitude lower than deduced from "standard" 
chemical models without depletion. This affects the dynamical evolution of prestellar cores by 
shortening the ambipolar diffusion time scale. 
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